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ABSTRACT 

'Evaluation  of  derivatives  and  Taylor  coefficients  of  functions  defined  by  computer 
programs  has  many  applications  in  scientific  computation.  The  process  of  automatic 
differentiation  of  such  functions  has  as  its  goal  the  production  of  machine  code  for 
the  evaluation  of  derivatives  and  Taylor  coefficients.  This  is  in  contrast  to  symbolic 

differentiation,  where  the  desired  output  is  a  more  or  less  pretty  formula,  and  numer- 

/ 

ical  differentiation,  which  is  inaccurate  and  unstable.  Another  distinction  between 
automatic  and  symbolic  differentiation  is  that  the  latter  usually  involves  considerable 
computational  overhead,  while  automatic  differentiation  can  be  carried  out  at  compile 
time  by  a  compiler  which  permits  user-defined  data  types  and  operators.  This  report 
shows  how  PASCAL-SC,  a  compiler  of  this  type,  can  be  used  to  generate  the  real  deriva¬ 
tive  types  GRADIENT,  HESSIAN,  TAYLOR,  and  the  corresponding  interval  types  IGRADIENT, 
IHESSIAN,  ITAYLOR.  Applications  of  these  types  to  solution  of  systems  of  nonlinear 
equations,  sensitivity  analysis,  constrained  and  unconstrained  optimization,  and  the 
solution  of  initial-value  problems  for  systems  of  ordinary  differential  equations  are 
indicated. v 
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SIGNIFICANCE  AND  EXPLANATION 


Many  problems  in  applied  mathematics  cam  be  solved  by  the  use  of  derivatives 
and  Taylor  series.  While  such  methods  appear  to  be  effective  in  theory,  their  appli' 
cation  in  scientific  computation  has  been  limited  until  recently  by  the  need  to  code 
the  required  derivatives  by  hand,  even  though  differentiation  is  a  well-understood 
mechanical  process.  Thus,  the  user  of  such  methods  has  had  the  choice  of  resorting 
to  numerical  differentiation,  am  inaccurate  and  unstable  process,  or  lately  to  the 
use  of  codes  for  symbolic  operations  which  include  differentiation.  The  latter  tend 
to  be  large,  unwieldy  systems,  intended  to  produce  formulas  and  perform  a  number  of 
symbolic  manipulations  in  addition  to  differentiation.  Ohis  report,  on  the  other 
hand,  describes  the  use  of  a  modern  compiler  which  permits  user-defined  data  types 
and  operators  in  order  to  perform  automatic  differentiation  of  computed  functions  in 
a  fast,  neat  way.  Here,  the  compiler  itself,  and  not  some  outside  system,  does  the 
differentiation  and  generation  of  Taylor  coefficients  desired,  and  produces  compact, 
efficient  machine  code  for  their  evaluation.  The  results  are  immediately  applicable 
to  problems  including  the  solution  of  systems  of  nonlinear  equations,  sensitivity 
analysis,  constrained  and  unconstrained  optimization,  the  solution  of  initial-value 


lies  with  MRC,  and  not  with  the  author  of  this  report. 
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1.  Automation  of  evaluation  and  differentiation  of  functions.  The  value. 

f(x)  of  a  function  of  n  variables,  x  ■  (x, ,...,x  )  is  obtained  on  a  computer  by 

I  n 

the  composition  f  -  f  of  o. . .of  of  a  finite  number  of  functions ,  each  usually  very 
simple,  such  as  an  arithmetic  operation  or  one  of  a  small  number  of  library  func¬ 
tions.  Using  the  fact  that  the  rule  for  evaluation  of  each  is  known,  it  is 
possible  for  a  computer  program  called  a  compiler  to  produce  machine  code  for  the 
evaluation  of  a  function  f  specified  in  some  form  similar  to  ordinary  mathematical 
notation,  and  including  the  possibility  of  different  cases  depending  on  x  and  inter¬ 
mediate  results.  A  key  feature  of  such  compilers  is  the  ability  to  perform  lofuruia 
translation,  that  is,  to  produce  machine  code  for  the  evaluation  of  an  exp4 MSton 
such  as 

(1.1)  F  <X*Y  +  SIN (X)  +  4 .0) * (3* (Y**2)  +  6); 
which  denotes  the  mathematical  function 

(1.2)  f(x,y)  -  (xy  +sinx  +  4)  (3y2  +  6). 

The  same  strategy  can  be  applied  for  the  machine  evaluation  of  derivatives  f ' (x) 
or  Taylor  coefficients  offatx  [6],  [9],  [10],  [11],  [12],  [13].  The  chain-rule 
of  calculus  applied  to  f  gives 

(1.3)  f'(x)  -  f*  (x(B_l))«...-f'  (xU))*f'  (x)  , 

x  m-l  m 

where  the  1  in  general  denotes  Frdchet  differentiation  and  the  •  matrix  multiplication 


(1.4) 


x(k)  -f<k)(x),  f(k) 


m-k+1 


►  • • 0 (  ,  k  =  1, . . . |B*1 • 
m 


Since  the  rules  for  calculation  of  derivatives  and  Taylor  coefficients  for  the  basic 
arithmetic  operations  are  perfectly  well  understood,  explicit  expressions  are  avail¬ 
able  for  the  functions  f|  [13] ,  and  thus  machine  code  can  be  generated  easily  for 
the  evaluation  of  derivatives  or  Taylor  coefficients  of  f  at  points  x  for  which  it 


is  differentiable  or  analytic,  respectively.  Most  existing  compilers,  however,  do 
not  present  this  opportunity  to  the  user.  The  reason  is  that  most  were  developed 
only  to  work  with  integers  and  floating-point  numbers  (i.e.,  the  types  INTEGER  and 
REAL)  originally,  and  are  quite  limited  in  the  kinds  of  data  and  operations  which 
can  be  handled. 

PASCAL-SC  [  3 ] ,  on  the  other  hand ,  was  developed  with  the  needs  of  scientific 


computation  in  mind,  and  thus  provides  complex  numbers,  intervals,  complex  intervals, 
and  vectors  and  matrices  over  these  types,  as  well  as  the  appropriate  operations,  in 
addition  to  integers,  and  floating-point  numbers,  vectors,  and  matrices  [17].  This 
vast  improvement  over  the  ordinary  type  of  compiler  is  achieved  by  allowing  the  in¬ 
troduction  of  user-defined  data  types  and  operators.  This  facility,  as  will  be  ex¬ 
plained  below,  permits  automatic  differentiation  and  generation  of  Taylor  coefficients, 
simply  by  introduction  of  the  appropriate  derivative  data  types  and  corresponding  op¬ 
erators  . 

2.  Derivative  data  types.  The  first  derivative  of  a  function  f  of  n  variables 


x  ■  (x1,...,xn)  is  its  gAadiznt  vtctoK 


(2.1) 


Vf(x)  -  ( 


3f  (x) 


3f  (x) 


similarly,  the  second  derivative  of  f  at  x  is  represented  by  its  He&Mdn  mntJux 


(2.2) 


Hf(x)  -  (  f-  f;fX)  )  , 


a  symmetric  n*n  matrix  [11].  Furthermore,  if  f  has  a  convergent  Taylor  series  ex¬ 
pansion  at  x  -  xQ,  then 


.  '•*  k  ¥  «  *  P 


(2.3) 


f(x)  -  I  <k)  <xc>  (*  -  xQ)k  +  VVX'V  ' 

k-0 

where,  of  course,  ^yf(0) (xQ)  ”  f(xQ)  I11! •  1°  the  scalar  case,  x  is  a  single  real 

variable ,  and  f(x)  is  approximated  by  the  TayloK  potynonUat 


n-1 


ViV^-V  ■  J0  ki 


V,(k) 


(xQ)  (x  - 


(2.4) 

which  in  turn  can  be  represented  by  the  n -vector 


V 


k 


(2.5)  Pn(xQ,t)  -  (f1#f2,...,fn) 

of  nofmaJLized  Taylox  coe^-c (UentA 

<2-61  fi  ■«^7Tf<l’1,(xo)ti’1'  1  '  x  -  V 

of  f.  Thus,  in  the  real  case,  real  vectors  and  matrices  are  required  for  the  rep¬ 
resentation  of  derivatives  and  Taylor  coefficients.  Similarly,  interval  vectors  and 
matrices  will  be  necessary  to  represent  the  same  quantities  for  interval-valued 
functions.  These  types  are  all  provided  in  PASCAL-SC  [17].  The  declaration  for 
vectors  in  the  real  case  (RVECTOR)  runs  as  follows: 

CONST  DIM  *  ; 

(2.7)  TYPE  DIMTYPE  =  1.. DIM; 

RVECTOR  »  ARRAY[DIMTYPE]OF  REAL; 


To  declare  a  real  matrix  (RMATRIX) ,  one  adds 

(2.8)  TYPE  RMATRIX  *  ARRAY[DIMTYPE]OF  RVECTOR; 

to  the  above.  These  declarations  will  be  basic  to  the  following  discussion. 

2.1.  Type  GRADIENT.  The  basic  derivative  type  treats  the  value  f(x) 
of  f  at  x  and  its  gradient  vector  Vf(x)  at  the  same  point  as  the  entity  (f (x) ,7f (x)) , 


that  is,  the  pair  consisting  of  the  real  value  of  f  at  x,  and  its  n-dimensional 
gradient  vector  of  partial  derivatives.  For  this  purpose,  the  RECORD  data  structure 
of  PASCAL-SC  is  ideal.  To  declare  this  type,  (2.7)  is  followed  by 


(2.9) 


TYPE  GRADIENT  -  RECORD  F:  REAL;  DF:  RVECTOR  END; 


so  that  if  V  is  a  variable  of  type  GRADIENT,  then  V.F  will  be  its  value,  and  V.DF 
its  gradient  vector.  Thus,  if  the  function  v  is  denoted  in  the  computer  program  as 
V,  then  V.F  ■  v(x)  and  V.DF  ■  Vv(x)  at  the  current  value  x  of  the  independent  var¬ 
iables  x, ,...,x  .  Thus,  V.DF [i]  -  3v(x)/3x. ,  i  -  l,...,n.  The  ith  independent  var- 
X  n  x 

iable  will  be  denoted  by  a  GRADIENT  variable,  say  XI,  such  that  XI. F  =  x^,  the 
current  value  of  x^,  and  XI. DF  is  the  ith  unit  vector  e^  =  (0,. .. ,0,1,0,. .. ,0) 
which  has  its  ith  component  equal  to  1,  and  the  others  equal  to  zero.  If  it  is 
desired  to  represent  a  constant  c  as  a  GRADIENT  variable  C,  then  C.F  *  c,  while 
C.DF  will  be  the  zero  vector  (0,...,0).  The  user  of  type  GRADIENT  is  free  to  name 
and  order  the  independent  variables  arbitrarily. 

2.2.  Type  HESSIAN.  This  type  can  be  considered  to  be  an  extension 
of  type  GRADIENT.  Here,  operations  are  performed  on  the  triple  (f (x) ,Vf (x) ,Hf (x)) 
as  the  basic  datum.  Once  again,  the  RECORD  structure  is  appropriate,  and  the  dec¬ 
laration  of  this  type  consists  of  (2.7)  followed  by  (2.8)  and 

(2.10)  TYPE  HESSIAN  *  RECORD  F:  REAL;  DF:  RVECTOR;  HF:  RMATRIX  END; 

with  now  V.F  -  v(x) ,  V.DF  =  Vv(x) ,  and  V.HF  »  Hv(x)  for  a  variable  V  of  type  HESSIAN 

2 

corresponding  to  the  function  v.  Here,  V.HF[i  ,  j]  *  3  v(x) /Sx^Sx^ ,  for  example. 

2.3.  Type  TAYLOR.  As  indicated  above,  the  quantity  to  be  computed  with 
in  this  case  is  the  vector  (2.5)  of  normalized  Taylor  coefficients  (2.6).  It  is 
important  in  most  applications  to  be  aware  of  the  &caZ&  ^acJton.  t  ■  x  -  x  ,  which 
will  be  of  type  REAL  in  the  scalar  case.  Thus,  the  RECORD  structure  will  also  be 
used  for  this  derivative  type,  which  is  declared  by  (2.7)  followed  by 

(2.11)  TYPE  TAYLOR  =  RECORD  T:  REAL;  TC:  RVECTOR  END; 

with  V.T  -  x  -  xQ  and  V.TC[i]  -  |v(1~1)  (xQ)  (x  -  x^*’1  for  i  -  l,...,n.  The 

value  of  the  Taylor  polynomial  (2.4)  can  be  obtained  very  simply  and  accurately  in 
PASCAL -SC  by  use  of  the  SUM  function: 


(2.12) 


P  :=  SUM(V.TC,0); 


where  P  is  of  type  SEAL.  In  (2.12),  the  sum  is  rounded  to  the  closest  floating¬ 
point  number;  upward  or  downward  rounding  can  be  achieved  by  replacing  the  0  by  +1 
or  -1,  respectively  [17]. 

2.4.  Interval  types  I  GRADIENT.  IHESSIAN,  ITAYLOR.  Interval  arith¬ 
metic,  including  operations  for  interval  vectors  and  matrices  and  standard  functions, 
is  another  convenient  feature  of  PASCAL-SC  [ 3 ] ,  [17] .  By  the  use  of  interval 
computation,  inclusions  of  the  real  values  of  expressions  such  as  (1.1)  can  be  ob¬ 
tained  [  9 ] ,  [10] .  The  same  applies  to  values  of  derivatives  and  Taylor  coefficients , 
so  the  derivative  types  introduced  above  can  also  be  defined  over  intervals.  The 
standard  declaration  of  an  interval  in  PASCAL-SC  is: 


(2.13) 


TYPE  INTERVAL  *  RECORD  INF.SUP:  REAL  END; 


and  interval  vectors  and  matrices  are  declared  by 


(2.14) 


TYPE  IVECTOR  »  ARRAY[OIMTYPE]OF  INTERVAL; 
IMATRIX  -  ARRAY[DIMTYPE]0F  IVECTOR; 


corresponding  to  (2. 7) -(2. 8).  The  INTERVAL  versions  of  the  REAL  derivative  types 
GRADIENT,  HESSIAN,  and  TAYLOR  are  thus  declared  by 

TYPE  I GRAD I ENT  »  RECORD  IF: INTERVAL;  IDF: IVECTOR  END; 

(2-15)  IHESSIAN  =  RECORD  IF: INTERVAL;  IDF:IVECTCr.;  IHF: IMATRIX  END; 

ITAYLOR  *  RECORD  IT: INTERVAL;  ITC: IVECTOR  END; 

respectively. 

3.  Derivative  operators.  In  order  for  expressions  to  be  evaluated  correctly 
when  derivative  data  types  appear,  the  arithmetic  operations  and  standard  functions 
have  to  be  defined  in  such  a  way  as  to  incorporate  the  appropriate  rules  for  differ¬ 
entiation  or  generation  of  Taylor  coefficients.  Furthermore,  as  in  (1.1),  provision 
for  variables  or  constants  of  type  INTEGER  or  REAL  must  be  made.  The  following  rule 


applies : 


VaAiablu  ok  expKU&ton&  o$  type,  real  ok  integer  mJU  be  tKeated  cu  constants 
<$04.  the  puKpo&e  o£  di^eKenttcution. 

In  order  to  simplify  the  following  discussion,  generic  variables  of  type  INTEGER 
will  be  denoted  by  K,  and  of  type  REAL  by  R,  RA,  RB.  The  derivative  types  introduced 
in  S2  will  be  indicated  by  G,  H,  T,  IG,  IH,  IT,  respectively,  and  the  general  deriva¬ 
tive  type  by  D.  Thus,  if  D  *  T,  type  TAYLOR  is  meant.  Generic  variables  of  type  D 
will  be  denoted  by  D,  DA,  DB  for  D  €  {G,H,T,IG,IH,IT}.  The  necessary  arithmetic 
operators  are  given  below;  operations  between  reals  and  types  IG,IH,IT  are  undefined. 

3.1.  Addition  operators. 

(3.1)  +D,  K+D,  D+K,  R+D,  D+R,  DA+DB. 

3.2.  Subtraction  operators. 

(3.2)  -D,  K-D,  D-K,  R-D,  D-R,  DA-DB. 

3.3.  Multiplication  operators. 

(3.3)  K*D ,  D*K,  R*D,  D*R,  DA*DB. 

3.4.  Division  operators. 

(3.4)  K/D,  D/K,  R/D,  D/R,  DA/DB. 

There  are  thus  22  operators  for  each  real  type  and  14  for  each  interval 
derivative  data  type.  For  example,  taking  D  ■  H,  the  source  code  for  the  operator 
symbolized  by  H+R  (addition  of  a  REAL  to  a  HESSIAN)  is: 

OPERATOR  +  (H:  HESSIAN ;R:  REAL)  RES:  HESSIAN; 

VAR  U:  HESSIAN; 

(3.5)  BEGIN  U.F:=H.F+R;U.DF:=H.DF;U.HF:=H.HF; 

RES:=U 

END; 

since  the  addition  of  a  constant  alters  only  the  value  of  a  variable,  and  not  its 
derivatives.  Rules  for  performing  the  above  operations  for  each  derivative  type 


are  given  in  [13] ,  for  example .  Source  code  for  each  of  the  real  data  types  G,H,T 
will  be  given  for  the  operation  DA*OB  in  the  next  section.  A  complete  set  of  source 
code  for  arithmetic  and  power  operations  and  standard  functions  can  be  found  in  [14] 
for  type  GRADIENT. 

3.5.  Power  operators.  The  provision  of  the  power  operator  **  to 
evaluate  xy  is  not  standard  in  PASCAL -SC.  In  order  to  implement  this  operator  for 
derivative  data  types,  it  is  necessary  to  define  it  for  some  combinations  of  the 
basic  REAL  and  INTERVAL  types.  Thus,  one  needs 

(3.6)  R**K,  RA**RB, 

for  type  REAL,  and,  if  I,IA,IB  denote  generic  INTERVAL  variables , 

(3.7)  I**K,  K**I,  IA**IB 

for  type  INTERVAL.  In  terms  of  these  basic  power  operations,  the  derivative  opera¬ 
tors  symbolized  by 

'3.8)  D**K,  K**D,  D**R,  R**D,  DA**DB 

can  be  defined.  Thus,  5  operators  for  each  of  the  derivative  types  are  needed,  in 
addition  to  the  2  operators  (3.6)  for  type  REAL,  and  the  3  operators  (3.7)  for  type 
INTERVAL,  as  well  as  D**K,  K**D,  DA**DB  for  D  6  (IG,IH,IT>. 

The  operator  R**K  can  be  implemented  by  repeated  squaring  [13] ,  [14] ,  or, 
better  yet,  by  the  accurate  algorithm  described  by  Btthm  [ 2 ].  RA**RB  can  be  cal¬ 
culated  using  R**K  for  the  integer  part  of  RB,  and  the  ordinary  exponential  func¬ 
tion  for  the  logarithm  of  the  base  multiplied  by  the  fractional  part  of  rb.  For 

Y 

interval  variables  X,Y,  X  is  defined  by 

(3.9)  xY  ■  [min{xy  |  xSx,y€y> ,max{xy  |  xex,y€y}], 

and  the  operator  **  is  defined  accordingly.  Since  the  inclusion 

2  2  2 
X  -  [min{x  I  a  s  x  £  b},max{x  |  a  £  x  s  b}]  r  .  X 


13.10) 


can  be  proper  for  intervals  X  =  [a,b] ,  the  standard  function  ISQR  for  interval 
arithmetic  [17]  is  used  to  compute  1**2,  and  I**K  by  repeated  squaring. 

The  complete  package  of  arithmetic  operations  for  derivative  data  types  thus 
consists  of  27  operators  for  each  real  type,  the  two  real  operators  (3.6),  the  three 
interval  operators  (3.7)  ,  and  the  17  arithmetic  operators  for  the  interval  types. 

The  priorities  of  the  operators  given  in  this  section,  from  highest  to  lowest, 
are:  1°.  Unary  addition  and  subtraction,  symbolized  by  +D  and  -D,  respectively? 
2°.  Multiplication,  division,  and  power:  *,/,**%  3°.  Binary  addition  and  sub¬ 
traction  ±.  Persons  familiar  with  compilers  in  which  **  has  higher  priority  than 
*,/  should  beware,  and  use  parentheses ,  as  in  (1.1). 

4.  Examples  of  multiplication  operators  in  PASCAL-SC.  In  order  to  give  ex¬ 
plicit  examples,  the  source  code  for  definition  of  the  operator  *  for  the  real 
derivative  types  GRADIENT,  HESSIAN,  and  TAYLOR  will  be  given  in  this  section,  both 
in  component  form  and  compact  form  using  the  vector  and  matrix  operators  available 
in  PASCAL-SC.  Complete  source  code  for  the  29  operators  required  for  type  GRADIENT 
is  given  in  [14] . 

4.1.  GA*GB  for  type  GRADIENT.  In  component  form,  the  source  code  for 
this  operator  is: 

OPERATOR  *  (GA.GB:  GRADIENT)  RES:  GRADIENT; 

VAR  I:  DIMTYPE;U:  GRADIENT; 

BEGIN  U.F:sGA.F*GB.F; 

(4.1) 

FOR  I:-l  TO  DIM  DO  U.DF[I]:=GA.F*GB.DF[I]+GB.F*GA.DF[I]; 
RES:»U 

END; 

however,  since  multiplication  of  vectors  by  real  numbers  and  addition  of  vectors 
are  operations  which  are  available  in  PASCAL-SC  [17]  ,  these  can  be  used  to  produce 
the  more  compact  source  code 


OPERATOR  *  (GA,GB:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT; 

BEGIN  U . F : =GA . F*GB . F ; 

(4.2)  U.DF:=GA.F*GB.DF+GB.F*GA.DF; 

RES:=U 

END; 

4.2.  HA*HB  for  type  HESSIAN.  If  XI  is  the  ith  independent  variable 
of  type  HESSIAN ,  then  XI. DF  will  be  the  ith  unit  vector,  as  for  type  GRADIENT,  and 
XI. HF  will  be  the  zero  matrix.  If  a  constant  C  is  declared  as  type  HESSIAN,  then 
both  the  vector  part  C.DF  and  the  matrix  part  C.HF  of  C  will  be  zero.  It  is,  of 
course,  more  economical  to  introduce  constants  as  type  INTEGER  or  REAL.  Multipli¬ 
cation  of  HESSIAN  variables  HA,HB  is  accomplished  by  the  operator 

OPERATOR  *  (HA,HB:  HESSIAN)  RES:  HESSIAN; 

VAR  I.J:  DIMTYPE;U:  HESSIAN; 

BEGIN 

U.F:=HA.F*HB.F; 

FOR  I:=l  TO  DIM  DO 

U.DF[I] :«HA.F*HB.DF[I]+HB.F*HA.DF[I] ; 

FOR  I :*1  TO  DIM  DO 

(4.3)  FOR  J:=I  TO  DIM  DO 

BEGIN 

U.HF[I  ,  J]:=HA.F*HB.HF[I  ,  J]+HA.DF[I]*HB.DF[J]+ 
HB.DF[I]*HA.DF[J]+HB.F*HA.HF[I  ,  J]; 

IF  I <  >  J  THEN  U.HF[J  ,  I]:=U.HF[I  ,  J] 

END; 

RES:=U 

END; 

in  component  form,  where  the  well-known  formulas  for  the  first  and  second  deriva¬ 
tive  of  products  have  been  used  [ 7 ] ,  [13] .  The  calculation  of  u.DF  here  is  taken 
directly  from  (4.1).  The  first  FOR  loop  could  be  eliminated  by  computing  U.DF [I] 
in  the  second.  As  before,  (4.3)  can  be  simplified  by  using  the  vector  and  matrix 
operations  of  PASCAL-SC,  augmented  by  the  oatZA.  product  VA**VB  of  vectors  defined 
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>y  the  operator 


OPERATOR  **  (VA.VB:  RVECTOR)  RES:  RMATRIX; 

VAR  I , J:  DIMTYPEjU:  RMATRIX; 

BEGIN 

4.4)  FOR  I.—  1  TO  DIM  DO 

FOR  J:=l  TO  DIM  DO  U[I  ,  J]:-VA[I]*VB[J]; 

RES:=U 

END; 

or  example.  With  this  operation,  (4.3)  becomes 

OPERATOR  *  (HA.HB:  HESSIAN)  RES:  HESSIAN; 

VAR  U:  HESSIAN; 

BEGIN 

„  ..  U.F:=HA.F*HB.F; 

4.5) 

U . DF : =HA . F*HB . DF+HB . F*HA . DF ; 

U.HF:=HA.F*HB,HF+HA.DF**HB.DF+HB.DF**HA.DF+HB.F*HA.HF; 

RES:=U 

END; 

hich  is  more  compact  but  less  efficient  than  (4.3) 

4.3.  TA*TB  for  type  TAYLOR.  Here,  the  well-know  formula  for  the  Taylor 
oefficients  of  the  product  (9],  [13]  are  used.  In  symbolic  form,  for  U:=TA*TB, 
ne  has  ( [13] ,  p.  41) 

J 

4.6)  U.TC  [J]  =  l  TA.TC[I]*TB.TC[J-I+1] ,  J«1,...,DIM. 

1=1 

very  important  advantage  of  the  use  of  PASCAL-SC  for  this  calculation  is  that 
4.6)  can  be  computed  as  the  scalar  product  of  two  vectors;  this  is  done  by  the 
tandard  function  SCALP  to  the  closest  floating-point  number,  or  can  be  rounded  up 
c  down  if  desired.  This  means,  for  example,  that  U.TC [1] ,... ,U.TC [DIM]  can  each 
e  computed  with  only  one  rounding  error  from  the  components  of  TA.TC  and  TB.TC, 
istead  of  having  increasing  rounding  error  as  J  increases.  This  helps  to  ameliorate 
le  of  the  possible  bugaboos  of  the  use  of  Taylor  series  methods  in  scientific 
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computation.  As  can  be  seen  in  [13] ,  a  number  of  other  formulas  for  the  recursive 

generation  of  Taylor  coefficients  take  the  form  of  scalar  products,  so  that  one  can 

take  advantage  of  the  accuracy  of  SCALP  in  their  computation.  Source  code  for  the 

multiplication  operator  *  for  TAYLOR  variables  is 

OPERATOR  *  (TA,TB:  TAYLOR)  RES:  TAYLOR; 

VAR  I ,J,K:  DIMTYPE;X,Y:  RVECT0R;U:  TAYLOR; 

BEGIN 

IF  TA.T  <>  TB.T  THEN 
BEGIN 

WRITELNC ERROR:  MULTIPLICATION  OF  TAYLOR  VARIABLES  WITH 
UNEQUAL  SCALE  FACTORS' );SVR(0) 

END; 

?  U.T:aTA.T;F0R  I:«l  TO  DIM  DO 

BEGIN  X[I]:*0;Y[I]:*0  END; 

FOR  I :=1  TO  DIM  DO 

BEGIN  X[I]:=TA.TC[I]; 

FOR  J:=1  TO  I  DO 
BEGIN 

K:*I-J+T;Y[J]:*TB.TC[K] 

END; 

U.TC[I]:=SCALP(X,Y,0) 

END; 

RES:=U 

END; 

where  an  attempt  to  multiply  TAYLOR  variables  with  unequal  scale  factors  results  in 
printing  an  error  message  and  return  of  control  to  the  operating  system.  The  in¬ 
itialization  of  the  vectors  X,Y  could  be  done  by  the  assignments  X : “VRNULL ,  Y : -VRNULL , 
from  a  parameterless  function  available  in  PASCAL-SC  to  produce  zero  vectors  [17]. 

The  present  implementation  of  type  TAYLOR  is  for  the  scalar  case,  with  a  single 
independent  real  variable,  say  x,  and  expansions  sure  performed  at  x  ■  xQ.  In  this 
case,  the  corresponding  independent  variable  X  of  type  TAYLOR  is  defined  as  follows: 

(4.8)  X.T  -  X  -  xQ,  X.TC [1]  -  xQ,  X.TC[2]  -  x  -  X  ,  X.TC[I]  -  0  for  I  >  2. 
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Thus,  the  assignment 


(4.9)  X.TC[1]:=X.TC[1]+X.TC[2];  or  X.TCp ]:=SUM(X.TC,0) ; 

will  generate  a  sequence  x^  -  xi_1  +  h,  h  =  x  -  xQ  of  points  of  expansion  of  the 
corresponding  TAYLOR  variables  dependent  on  X. 

Multiplication  of  the  interval  derivative  types  I GRADIENT,  IHESS1AN,  AND  ITAYLOR 
follows  the  pattern  given  above  for  the  corresponding  real  types. 

5.  Standard  derivative  functions.  The  rules  for  differentiation  of  the  standard 
functions  ABS,  SQR,  SQRT,  EXP,  Hi,  SIN,  COS,  ARCTAN  are  well  understood,  as  is  the 
generation  of  Taylor  coefficients  for  these  functions.  Thus,  there  is  no  conceptual 
(or  practical)  difficulty  in  providing  the  standard  functions  DABS,  DSQRT,  DEXP, 

DIN,  DSIN,  DCOS ,  DARCTAN  for  the  derivative  types  D  6  {G,H,T,IG,IH,IT) ,  and  the 
function  DSQR  for  the  interval  derivative  types  D  6  {1G,IH,IT}. 

Naming  functions  in  this  way  mans  that  in  an  expression  for  the  function  (1.2) , 
Taylor  coefficients  are  obtained  automatically  if  one  declares  F,X,Y  to  be  of  type 
TAYLOR,  and  writes  the  assignment 

(5.1)  F  :=  (X*Y  +  TSIN(X)  +  4.0) * (3* (Y**2)  +  6); 

here,  the  minor  nuisance  of  writing  TSZN  instead  of  SIN  is  offset  by  the  fact  that 
the  former  makes  it  clear  that  the  expression  in  (5.1)  is  of  type  TAYLOR,  and  that 
arguments  of  type  TAYLOR  are  expected. 

A  complete  set  of  source  code  for  the  GRADIENT  (D  *  G)  versions  of  the  standard 
functions  listed  above  is  given  in  [14].  For  example,  for  the  gradient  cosine, 

FUNCTION  GC0S(G:  GRADIENT):  GRADIENT; 

VAR  M:  REAL;U:  GRADIENT; 

(5.2)  BEGIN 

M:*-SIN(G.F) ;U. F:*C0S(G.F) ;U.DF:sM*G.DF;GCOS:»U 
END; 

in  vector  form.  With  the  arithmetic  operations  and  the  standard  functions  given,  it 
is  very  easy  for  the  user  to  introduce  others.  If  the  gradient  secant  is  needed,  one 
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can  simply  use  the  code 


FUNCTION  GSEC(G:  GRADIENT):  GRADIENT; 

(5.3) 

BEGIN  GSEC:=1/GC0S(G)  END; 

for  this  purpose,  and  so  on. 

6.  Applications  of  derivative  types  in  scientific  computation.  Since  the 
operations  of  differentiation  and  series  expansion  are  ubiquitous  in  mathematical 
modeling,  possible  applications  of  derivative  data  types  in  scientific  computation 
would  fill  a  vast  catalog,  which  would  also  probably  never  be  completed.  Here, 
only  a  few  applications  which  have  been  implemented  in  software  in  one  form  or 
another  will  be  mentioned.  In  each  case,  the  PASCAL-SC  formulation  of  these  types  and 
their  corresponding  functions  and  operations  leads  to  a  drastic  simplification  of 
previous  efforts. 

6.1.  Applications  of  type  GRADIENT  and  IGRADIENT.  A  simple  application 
of  type  GRADIENT  is  to  anaZyiZi .  If  F  is  of  type  GRADIENT,  then  F.DF  [i] 

•  3f/3x.  is  the  rate  of  change  of  the  function  f  symbolized  by  F  with  respect  to  the  ith 
independent  variable  .  Similarly,  F.DF [i]/F.F  is  the  XZJtative.  rate  of  change  of  f 
with  respect  to  x^ .  Furthermore,  the  gradient  vector  F.DF  -  Vf (x)  gives  the  direction 
of  the  fastest  increase  of  f  at  the  current  values  of  the  independent  variables,  a 
fact  which  is  useful  in  optimization  and  other  applications . 

Given  several  GRADIENT  variables  F,G,H,  their  Jacobean  meutrux  J  will  have  rows 
J[l]  -  F.DF,  J[2]  =  G.DF,  J[3]  -  H.DF.  Knowing  the  values  F.F,  G.F,  and  H.F  of  these 
variables  as  well  as  their  Jacobian  matrix  makes  it  very  easy  to  code  the  solution 
of  systems  of  equations  by  Newton's  method  in  terms  of  GRADIENT  variables  [ 7 ] ,  [11] , 
[14].  Furthermore,  the  type  IGRADIENT  can  be  used  to  obtain  Lipschitz  constants  for 
functions  of  several  variables  [14] . 

In  interval  analysis  [ 9 ] ,  [10] ,  it  is  well-known  that  the  mean-value  form 
(6.1)  F (X)  -  f(x)  +  F'(X)(X  -  x) 

can  give  an  accurate  interval  inclusion  of  a  real  function  f  on  a  small  interval 
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X  containing  x.  If  f  is  a  function  of  several  variables,  an  interval  inclusion 
F' (X)  of  the  Jacobian  matrix  f 1 (x)  is  needed.  This  can  be  obtained  automatically 


if  f  is  coded  n  terms  of  its  components  f^  as  variables  FI  of  type  I GRADIENT. 

By  use  of  PASCAL-SC  and  the  types  GRADIENT  and  I GRADIENT ,  the  program  NEWTON 
[  7 ]  for  the  solution  of  nonlinear  systems  of  equations  can  be  reduced  from  over 
3,400  lines  in  FORTRAN  to  a  few  dozen. 

6.2.  Applications  of  type  HESSIAN  and  IHESSIAN.  As  one  might  expect, 
these  types  appear  to  be  extremely  useful  in  optimization.  In  unconstrained  opti¬ 
mization,  the  gradient  vector  Vf(x)  of  a  function  f(x)  =  f(x. ,...,x  )  will  vanish 

I  n 

at  a  maximum  or  minimum  value  of  f ,  so  one  wants  to  solve  the  system  of  equations 
Vf(x)  ■  0,  the  Jacobian  matrix  of  which  is  Hf(x).  If  f  is  coded  as  a  HESSIAN  var¬ 
iable  F,  that  is,  by  an  assignment  of  the  form 

(6.2)  F  :=  F(X1,...,XN) ; 

then  the  components  of  the  gradient  vector  are  given  automatically  by  F.DF,  and 
the  Hessian  matrix  of  f  by  F.HF,  the  latter  being  useful  in  solution  procedures 
based  on  Newton's  method  [14].  Suppose  that  the  optimization  problem  is  con&tJWAJnzd 
by  conditions  g(xir...,xn)  =  0,  h(x^,...,xn>  ■  0,  where  g,h  are  coded  in  symbolic 
form  by  the  assignments 

(6.3)  G  :=  G(X1,...XN);  H  H (Xl , . . . ,XN) ; 

as  variables  of  type  HESSIAN.  Then,  increasing  the  dimension  of  the  problem  to  n+2 
and  introducing  two  new  independent  variables  Ll,  L2  (the  Lag/umge.  nuttipLLeJU)  as 
of  type  HESSIAN,  one  makes  the  assignment 

(6.4)  W  :»  F  +  Ll *G  +  L2*H; 

and  the  solution  of  the  unconstrained  problem  for  W  is  then  a  solution  of  the  con¬ 
strained  problem  for  F  [14].  Once  again,  the  Hessian  matrix  W.HF  of  W  is  the  Jacobian 
matrix  of  the  system  W.DF  =  0,  which  gives  a  necessary  condition  for  an  extreme  value 
of  W,  and  can  be  used  in  solution  procedures  based  on  Newton's  method. 


In  the  program  NEWTON  [ 7  ] ,  interval  Hessian  matrices  were  used  to  obtain 


Lipschitz  constants  for  Jacobian  matrices  and  ultimately  rigorous  error  bounds  for 
approximate  solutions  of  nonlinear  systems  of  equations.  These  matrices  can  be 
obtained  automatically  if  the  type  IHESSIAN  is  used  in  the  computation  [14] . 

6.3.  Applications  of  type  TAYLOR  and  ITAYLOR.  Taylor  series  and 
the  approximation  of  functions  by  Taylor  polynomials  appear  in  many  places  in 
scientific  computation.  One  of  the  most  successful  areas  of  application  of  the 
idea  of  recursive  generation  of  Taylor  coefficients  is  the  numerical  solution  of 
the  initial-value  problem  for  ordinary  differential  equations  and  systems  of  such 
equations  [  1  ] ,  [ 4  ] ,  [  9  ]  ,  [10] .  Type  TAYLOR  is  ideally  suited  for  production 
of  software  for  this  purpose.  For  example,  suppose  a  Taylor  polynomial  approximation 
is  sought  for  the  solution  of  the  equation 

(6.5)  y'  =  (xy  +  sinx  +  4)  (3y2  +  6),  y(xQ)  *  yQ» 

say  on  the  interval  xQ  £  x  <_  x^.  Part  of  the  PASCAL-SC  code  for  this  purpose, 
with  X,  Y,  YPRIME  of  type  TAYLOR,  would  run  as  follows  for  a  REAL  step- length  of 
H: 


|  BEGIN  X.T:-HjY.Ts«H;X.TC»-VRNULLjX.TC[l] :-X0jX.TC [2] :-H» 

*  WHILE  X.TCI1]  <=  XT  DO 

J  BEGIN  Y.TC:“VRNULL ? Y.TC [1] :-Yp; 

;  FOR  Ij-2  TO  DIM  DO 

j  BEGIN 

YPRIME s* (X*Y+TSIN (X)  +4) * (3* (Y**2) +6) ;  J;«I+1; 

(6.0) 

■  Y.TCtI] :-YPRIME.TC  [J] *H/J 

END; 

X.TC  [1]  :*X.TC [  1]  +H;  YJ8  :*SUM (Y .TC  ,0)  ; 

WRITELN (X.TC [1] ,Y 0) 

END; 

END; 

I 

The  same  kind  of  coding  can  be  expanded  for  systems  of  equations. 

In  the  prograun  DIFEQ  ( 4  ] ,  interval  Taylor  coefficients  and  polynomials  are 
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used  for  rigorous  error  and  step-size  control  in  the  computation  analogous  to 
(6.6)  for  the  recursive  generation  of  Taylor  coefficients  of  y  from  those  for  y'. 

In  the  program  INTE  [ 5  ] ,  interval  Taylor  coefficients  are  used  to  obtain  guaranteed 
bounds  for  the  error  terms  in  various  rules  of  numerical  integration.  In  both  of 
these  cases,  the  use  of  the  type  ITAYLOR  would  result  in  the  reduction  of  programs 
with  thousands  of  lines  in  FORTRAN  and  assembly  language  to  PASCAL-SC  programs  of 
a  hundred  lines  or  so. 

A  challenge  in  the  extension  of  type  TAYLOR  to  the  vector  case  is  that  deriva¬ 
tives  are  then  multilinear  operators  [11] ,  so  that  both  storage  requirements  and 
operation  times  increase  drastically.  However,  some  of  the  present  and  projected 
"supercomputers"  may  well  be  suitable  for  Taylor  series  methods  to  be  applied  to 
partial  differential  equations  along  the  lines  indicated  here,  given  compilers  with 
the  capabilities  of  PASCAL-SC  which  can  also  exploit  any  parallelism  available  in 
the  hardware . 
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puter  programs  has  many  applications  in  scientific  computation.  The  process  of 
automatic  differentiation  of  such  functions  has  as  its  goal  the  production  of 
machine  code  for  the  evaluation  of  derivatives  and  Taylor  coefficients.  This 
is  in  contrast  to  symbolic  differentiation,  where  the  desired  output  is  a  more 
or  less  pretty  formula,  and  numerical  differentiation,  which  is  inaccurate  and 
unstable.  Another  distinction  between  autosiatic  and  symbolic  differentiation 
is  that  the  latter  usually  involves  considerable  computational  overhead,  while 
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ABSTRACT  (continued) 


automatic  differentiation  can  be  carried  out  at  compile  time  by  a  compiler  which 
permits  user-defined  data  types  and  operators.  This  report  shows  how  PASCAL-SC, 
a  compiler  of  this  type,  can  be  used  to  generate  the  real  derivative  types  GRADIENT, 
HESSIAN,  TAYLOR,  and  the  corresponding  interval  types  I GRADIENT,  IHESSIAN,  ITAYLOR. 
Applications  of  these  types  to  solution  of  systems  of  nonlinear  equations,  sensi¬ 
tivity  analysis,  constrained  and  unconstrained  optimization,  and  the  solution  of 
initial -value  problems  for  systems  of  ordinary  differential  equations  are  indicated. 


